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Determining  Conductivity  by  Boundary  Measurements: 
Some  Numerical  Results 


In-Ja  B.  Lee| 


Abstract.  We  consider  the  problem  of  determining  an  unknown  conductivity  by  boundary  measure¬ 
ments.  First  we  briefly  review  various  known  results  about  uniqueness  and  continuous  dependence.  The 
main  emphasis  of  this  paper  is  on  the  discussion  of  some  numerical  results  for  a  variational  algorithm  to 
recover  the  conductivity. 


1.  Introduction.  We  study  the  following  inverse  problem:  can  one  determine  an 
unknown  conductivity  7(1)  inside  a  body  ft  by  means  of  measurements  of  potential  and 
flux  on  the  boundary?  This  problem  was  first  raised  by  Calderon  and  has  since  been 
studied  by  several  people. 

In  the  following  ft  will  be  a  bounded  C°°  domain  in  Rn,  n  >  2,  with  boundary  I\  The 
unknown  conductivity  7(2)  is  in  L°°(ft)  and  satisfies  0  <  70  <  7(2)  for  all  x  €  ft.  Consider 
the  operator 

I7(u)  =  V.(7(x)Vu) 

acting  on  if1  (ft). 

We  define  the  Dirichlet  to  Neumann-data  map 

A7  :  H1/2(T)  ->  iJ~1/2(r)  by 
du 

=  75T’ 

where  u  solves  the  following  boundary  value  problem 

L-/U  =0  in  ft, 

(1.1) 

u  =  d  on  T. 

The  inverse  problem  now  is  to  determine  7  given  knowledge  of  A7(<£).  We  say  that 
7  is  identifiable  if  the  map  7  1— >  A7  is  injective.  In  section  2  we  briefly  review  known 
results  concerning  the  identifiability  of  7.  Section  3  concerns  the  question  of  continuous 
dependence.  In  the  last  three  sections  of  this  paper,  we  study  a  particular  algorithm  to 
determine  7;  the  emphasis  is  on  the  presentation  of  various  numerical  experiments. 

tDepartmem  of  Mathematics,  University  of  Maryland,  College  Park,  Maryland  20742.  Supported  by 
NSF  grant  DMS-8G01490  and  ONR  contract  N00014-85-K-0169. 


2.  Identiflability  Results.  For  many  purposes  it  is  more  convenient  to  work  with 
the  energies 

Qi(<t>)  =  f  7  |Vu|2  dx, 

Jn 

instead  of  the  operator  A7.  It  is  easy  to  see  that  knowing  Qy(<i>),  for  any  <j)  6  is 

equivalent  to  knowing  A7. 

Kohn  and  Vogelius  [5]  proved  that  7  is  uniquely  determined  by  Qy(4>)  if  7  is  analytic. 

THEOREM (KOHN-VOGELIUs).  Let  7 i  =  1,2,  be  in  T°°(fl)  with  a  positive  lower 
bound.  Let  xo  €  T  and  let  B  be  a  neighborhood  of  xo  relative  to  Q.  Suppose 

7  ieC°°(B),  for  i  =  1,2 


and 


Qi\{4>)  =  whenever  <f>  6  H1^2(T)  with  supp  <f>  c  R  n  r. 


Then 


Dkh(x0)  =  Dk  12(20),  for  all  k  =  (klt  •  •  •,  &n)  >  0, 

where  Dk  denotes  the  derivative  (d/dx1)kl  •  •  •  (3/3x„)fcn. 


The  main  trick  in  the  proof  of  this  theorem  is  to  take  boundary  data  d>  that  are  highly 
oscillatory,  with  vanishing  moments  and  supported  in  a  small  neighborhood  of  x0  €  T.  For 
such  if),  the  solution  of  (1.1)  will  decay  rapidly  away  from  x0;  these  <f>  are  therefore  well 
suited  for  studying  7  at  x0.  The  following  three  lemmas  are  the  main  ingredients  of  the 
proof  of  the  theorem. 


LEMMA  1.  Let  M  be  any  positive  integer,  and  let  z  be  a  point  on  T.  There  exists  a 
sequence  {<^/v})v=i  C  C°°(r)  such  that 


II^A||i/2+<,r  <  Ct  N*,  t  >  —M 


ll<Mli/2,r  =  1> 

supp <f>jv  i  {z}  as  N  — *  00. 


:  Accession  For 


Fiv  A/0  and  let  u n  denote  the  solution  to  the  boundary  value  problem 


A<3cl 


n 


LEMMA  2.  Let  £1'  C  fi  with  dist  (z,£l')  >  0,  and  assume  that  7  is  in  C°°  in  a 
neighborhood  of  z.  Then 

||uN||i,ft'  <  CN~m ° 

for  all  N  >  1.  The  constant  C  depends  on  7,  Q',  z,  and  Mo,  but  not  on  N. 

Define 

p(x)  =  dist  (x,  T)  for  any  x  6  Cl. 

LEMMA  3.  Assume  that  7  is  C°°  in  a  neighborhood  of  z,  and  let  U  denote  any  neigh¬ 
borhood  of  z.  Given  l  >  0  and  e  >  0,  there  exists  a  constant  Cii(  >  0  such  that 

[  p(x)l\VuN\2dx>C,,(N~{n+()l 

Ju 

for  sufficiently  large  N . 

Sketch  of  the  proof  of  Theorem:  the  proof  proceeds  by  contradiction  using  the  three 
lemmas.  If  the  theorem  were  not  true,  we  may  suppose  that 

7i(*)-72  (x)>Cp(x)‘,  xEU, 

for  some  constant  C  >  0  and  some  neighborhood  U  of  a  point  2  6  T,  arbitrarily  close  to 
x0.  Choose  Mo  >  \nl.  Let  u'N,  i  =  1,2,  solve 

L^ulN-fi ,  u‘N|r  =  <AA- 


Then,  for  sufficiently  large  N , 


dx  > 

/  72|Vu^ 

|2  dx  + 

C  /  p(x)‘\Vu'N\7 

dx 

Ju 

Ju 

> 

[  72|Vu)v| 

|2  dx  - 

/  72|Vu)v|2dx 

+  c^,eiv-(n+e), 

Jn 

Jn\u 

> 

f  72|Vu2„| 

\2dx  - 

CN~2Mo  +  C,,tN' 

-(«+«)* 

Jn 

j 

[  7i|Vu^|2 

’n 

dx  > 

f  72|Vu^|2dx, 

Jn 

Since  2A/0  >  nl ,  we  get 


for  all  sufficiently  large  N,  and  this  is  a  contradiction.  [] 

Kohn  and  Vogelius  [7]  have  extended  their  study  to  piecewise  analytic  7  showing  that  it 
is  possible  to  determine  an  unknown  conductivity  interior  to  a  body  Cl  from  measurements 
of  the  potential  and  the  flux  on  the  boundary. 

Sylvester  and  Uhlmann  [10]  have  subsequently  shown  the  uniqueness  result  in  tho  mom 
difficult,  case  7  £  Cl}. 


THEOREM(Sylvester-Uhlmann).  Let  Q  C  R"  (n  >  3)  be  a  domain  with  smooth 
boundary.  Suppose  70  and  71  are  in  C°°(fi),  70, 71  >  0  in  Q  and 

QyoW  =  QyAV  G  H1/2(F) 


then 

7o  =  7i  in  fi- 

This  theorem  has  been  established  only  for  space  dimension  greater  than  or  equal  to 
three.  It  is  not  known  whether  the  same  result  holds  for  space  dimension  two.  But  in  that 
case  Sylvester  and  Uhlmann  [9]  have  obtained  a  local  uniqueness  result  for  coefficients  that 
are  nearly  constant  in  C°°(fl). 

So  far,  we  only  have  considered  scalar  coefficients  7  corresponding  to  isotropic  conduc¬ 
tors.  It  is  also  of  practical  importance  to  study  symmetric  matrix  valued  coefficients 

7 ij  =  7 ji  G 

A|£|2  ^  (l(x)£,0  <  tI£|2  for  all  £  G  ft  and  £  G  R". 

A 

In  this  case  we  can  not  expect  to  recover  the  full  matrix  {7^}.  If  two  matrix  conductivities 
7  and  6  are  related  by  change  of  variables 

1  x — '  d$i 

ht  cni))  = 

».j  *  3 


($  is  a  diffeomorphism  of  Q),  then  u  solves  (1.1)  with  coefficient  7  if  and  only  if  u *  =  1/0$  1 
solves  the  similar  problem  with  coefficient  6.  If,  in  addition, 


$(1)  =  x  for  x  G  T, 

then  u*  =  u  on  T,  and 

7VU  •  v  =  iSVu’*’  ■  v  on  r. 

It  therefore  follows  that  7  and  6  have  identical  boundary  measurements.  Concerning 
uniqueness  Kohn  and  Vogelius  [6]  have  proved  that  if  (n  —  1)  eigenvalues  and  eigenvectors 
of  7  are  known  then  the  last  eigenvalue  can  be  distinguished  by  boundary  measurements. 

3.  Continuous  Dependence  Results.  It  is  very  important  to  study  the  continuity 
properties  of  the  mapping  A7  — ►  7.  We  define  the  operator  norm 


l|A7i|,/2,-l/2  =  SUp 

07*0 


ha7wh  h-'/2(  n 


Sylvester  and  Hhlmann  [11]  have  recently  shown 


THEOREM (SYLVESTER-UHLM ANN).  Suppose  that  71  and  72  are  measurable  and 
satisfy,  for  some  A, 

0  <  y  <  7^  <  A  for  i  —  1,2 

A 

then 

liA7i  ~  A-y2  II1/2.-1/2  <  Ci  |(7i  -72|U~(n)- 

If  in  addition,  71  and  72  are  continuous,  then 

hi  ~72|U°°(D  <  C2||ATi  -  A7j||]/2, -1/2- 
Let  Bi  —  A7i  —  7jAi-  If,  7i  and  72  are  Lipschitz  continuous,  and  for  some  (3, 

|V7i|  <  ft  for  i  =  1,2, 


then 


11-Si  ~  -^2 1|  1/2,— 1/2  <  C3II71  —  72 1! (C2) 


and 

li 7i  ~  72||vYi.~(r)  <  C4(||Bi  -  B2||i/2i-i/2  +  ||A71  -  A72 1|  1/2,— 1/2 ) 

The  constants  C,  depend  only  on  Q,  A  and  3. 

This  result  also  can  be  proven  by  a  simple  modification  of  the  method  of  Kohn  and  Vo- 
gelius.  See  for  instance  Alessandrini  [1].  Alessandrini,  furthermore,  obtained  the  following 
result  under  the  assumption  that  7 ,,  i  =  1,2,  satisfies  a-priori  inequalities 


Y  <  7 i(x),  for  all  x  € 

A 

ll7i II *  =  1,2 


for  some  s  and  A,  s  >  n/2. 

THEOREM  (Alessandrini).  With  conditions  as  stated  above, 


Il7i  -72|U«=(n)  <  w(||A7l  -  A-y2 )|i/2,— 1/2), 
and  the  function  u>  is  such  that 


uj(t)  <  C  \log  t\  s,  Vt,  0  <  t  <  1/e, 
for  some  C,  6,  0  <  6  <  1,  depending  only  on  A,  n  and  s. 


This  result  is  somewhat  disappointing  since  it  predicts  such  a  weak  form  of  continuous 
dependence.  In  its  generality  the  result  may,  however,  very  well  be  the  best  possible. 


a 

■& 


t 

t 


( 

X 

X 

A 

A 


a 


It  becomes  quite  natural  to  analyze  continuous  dependence  for  classes  of  conductivities 
with  more  restricted  spatial  dependence.  For  practical  reasons,  it  is  important  to  seek 
continuous  dependence  estimates  in  terms  of  finitely  many  sets  of  boundary  measurements. 

Friedman  and  Vogelius  [4]  studied  conductivities  that  correspond  to  a  finite  number 
of  small  inhomogeneities,  of  extreme  conductivity,  imbeded  in  an  n-dimensional  reference 
medium.  The  reference  conductivity  7  is  assumed  to  be  C2+f3 .  Each  inhomogeneity  has 
the  form  z k  +  e  pk  B,  where  B  is  some  bounded  domain  in  Rn,  n  >  2,  with  0  6  B  and  dB 
of  type  C2+P.  The  points  {17,  •  •  •,  Zk}  belong  to  ft  and  satisfy 

| Zk  —  zj\  >  do  >  0,  Vj  ^  k  and 
dist  (^t,r)  >  do,  Vfc. 

Here  e  is  assumed  to  be  small  enough  that  the  sets  {2*  4-  epkB}  are  disjoint  and  that  their 
distance  to  Rn\ft  is  larger  than  ft  is  assumed  to  be  bounded  with  boundary  T  £  C2+/?. 
Let 

K 

=  (J(~*  +  ePkB) 

k=l 

denote  the  total  collection  of  inhomogeneities.  If  they  all  ha\'e  infinite  conductivity,  then 
the  voltage  potential  u(,  given  the  boundary  current  ip,  solves 

min  {-  [  7|Vu|2dx—  j  1 puds}. 

u£H'  (!J),Vu=0  in  2  Jq  Jr 

For  this  problem  to  have  a  unique  solution  we  require  that 

J  0  ds  =  0  and  J  u(  ds  =  0. 

Consider  two  arbitrary  collections  of  inhomogeneities 


=  (J(z*  4-  epkB)  and  w'  =  |J(4  +  tpkB), 


and  denote  by  ut  and  u[  the  corresponding  voltage  potentials  with  fixed  boundary  current 
»/>.  Let  U  denote  the  solution  to 


It  is  important  that  U  satisfy  the  non-degeneracy  condition  of  VU  ^  0. 


min 


1 


THEOREM (Friedman-Vogelius).  Let  T j  be  nonempty  open  subset  of  T.  There 
exist  constants  0  <  to,  8$  and  C  and  a  function  77(e),  linage)  =  0,  such  that  if  e  <  e  0  and 

e— *0 

e-"IK  -  u'||Loo(ri)  <  60,  then 

(1)  I\  =  K' ,  and,  after  appropriate  reordering, 

(2)  \zk  ~  41  +  I  Pk  ~  P'k  I  <  C  e_n|K  -  u'||L«.(ri)  +  77(e)  (1  <  k  <  K). 

The  constants  (q  ,  <50  and  C  and  the  function  77  are  independent  of  the  two  sets  of  inhorno- 
geneities. 

A  similar  result  can  be  obtained  for  non-conducting  inhomogeneities  as  well  as  for  the 
case  where  the  voltage  potential  xp  is  fixed  and  the  current  flux  xp  is  measured. 

4.  Numerical  Algorithm.  To  build  an  algorithm  to  reconstruct  the  conductivity 
we  consider  the  functional 


F(7,  n,cr)=  f  |7J/2Vu  —  7  1/2cr|2  dx. 

J  n 


It  is  clear  that  F{-),u,o)  is  non-negative  and  that  F(^,u,o)  =  0  if  and  only  if  7 !/2Vu  = 
7-1/2<7,  i.e.,  if  and  only  if  7V11  =  a. 

Let  5  be  the  set  of  (u,o)  characterized  by 

S  =  {(u,  a)  |  u  =  <p  on  T,  V  -cr  =  0  in  ft  cr  ■  u  =  xp  on  T}. 

Then  F  =  0  on  5  if  and  only  if  7,  u,  and  <7  satisfy 

V  ■  (7 Vu)  =  0  in  ft 
u  =  <p  on  T 

(4.2)  gu 

1-^  =  4’  on  r 

<T  =  7  Vu. 

In  that  case,  7  is  a  conductivity  consistent  with  the  measurements  <p  and  xp. 

For  V  •  <7  =  0  in  ft,  u  =  <p  and  0  •  v  =  xp  on  T,  a  simple  computation  gives 

F(y,u,a)  =  I  7I Vu\2  dx  +  I  7_1  \<r\2  dx  —  2  j  <pxl>ds. 

Ju  Jq  Jr 

If  more  than  one  set  of  boundary  measurements  {</>,,  V’t}1/£1  are  given,  then  the  corre¬ 
sponding  functional  would  be 

Kf 

(4.3)  F(7,  {or,}^!)  =  V  /  |71/2Vu,-7  1/2<7,|2dr. 

,=1 


S3 


« 


$ 


The  set  S  would  be  denned  by 

5  =  {(it,-, crj)fL1 1  u,  =  <t>i  on  T,  V  •  cri  =  0  in  ft,  at  ■  u  =  xpi  on  T}. 

The  reconstruction  algorithm  we  consider  seeks  to  minimize  the  functional  F  over  5. 
Since 

bC  K  K 

F(7,  {ui}^i,{CT,}^i)  =  ^  f  7|Vu,|2  dx  4-  ^  f  7_1|crt|2dx  -  2^  f  <f>i  V,  ds, 

.=i  Ja  i=!  Jr 

the  minimization  of  F(7,{ui}-L1,{(7i}fl1),  for  fixed  7,  is  equivalent  to 

min  /  7|Vu,'|2  dx  subject  to  tq|r  =  4>i 

Jo. 

and 

min  /  7 — 1  j | 2  dx  subject  to  V  •  a,  =  0,  <7,  •  v  =  ipi, 

Ju 

for  all  1  <  i  <  K. 

The  above  minimization  is  equivalent  to  finding  the  solutions  of 

V  •  (7 Vui)  =  0  in  Q 

(4.4) 

it,  =  <f>i  on  r 

and 

V  •  (7 Viq)  =  0  in  ft 

(4.5)  7^1  =  0,  on  T 

dv 

with  <7;  =  7  Vu,-, 

for  all  1  <  i  <  A'. 

For  fixed  }  and  the  minimization  over  7  £  can  be  done  pointwise 

on  the  integrand 

^(7|Vu,|2 +7_1|(7i|2). 

i=i 

The  result  is 

7  —  (  1  -  — — ,  provided  of  course  this  is  in 

Ei= 1  lVu>l2 

Alternating  Direction  Algorithm:  [8],  [12],  [13]. 

1.  Pick  an  initial  guess  7. 

2.  For  fixed  7, 

2. a  Solve:  V  •  (7 Vut)  =  0  in  ft,  u,  =  <^>,  on  T,  1  <  i  <  K. 

2.b  Solve:  V-(jVvt )  =  0  in  ft,  7  dvjdv  =  on  T  and  set  a ,  =  7 Vi>,,  1  <  /  <  A  . 


>  ^  A  « V-V  r.’  \  %'  *»’  •.*  O'.  *  v  V  %*  V  %’ 


“,l  >*•***>»'  ■  ■"**  w'"*1  *  *  (,  ^  V  w'v 


3.  Update  7  by  minimizing  F( 7,  {u.J/ij,  {cr,},^.,)  for  fixed  {u, } and  {crJ-L,. 

4.  If  the  reconstruction  is  satisfactory,  exit1  ;  else  go  to  step  2. 

This  algorithm  concerns  scalar  7  only.  To  build  a  similar  algorithm  to  recover  a  matrix 
conductivity,  the  only  change  we  have  to  make  concerns  the  update  of  7.  We  still  minimize 


■^(7,  {^.},=i)  =  XJ[^{^Vu*')TVu'  +  (7-1  (Ti)T(Ti}dx  -  J  0,0,  ds] 


for  fixed  and  {<r,}  ■'L,  but  the  formulas  are  slightly  more  complicated: 

3. a  Compute  L  =  (Vui)7’. 

3.b  Compute  M  =  (cr«)T- 

3.c  Compute  eigenvalues  of  ML,  si  and  s 2  and  corresponding  eigenvectors  ej,  and  e2, 
and  set  A,  =  \/si,  for  i  =  1,2. 

3.d  Set  E  =  [e  1,  e2]  and  let  A  be  the  diagonal  matrix  with  entries  Aj  and  A2. 

3.e  Compute  X  =  Eh.E~l . 

3.f  Set  7  =  XL-1 

This  describes  the  alogrithm  to  recover  a  scalar  as  well  as  a  matrix  conductivity. 

We  have  tested  this  algorithm  on  several  test  problems  and  obtained  somewhat  promis¬ 
ing  results.  In  subsequent  sections,  we  present  the  numerical  results  obtained  from  this 
algorithm. 

A.  Wexler  et  al.  [12],  [13]  originally  proposed  an  algorithm,  which  is  similar  to  the 
algorithm  presented  here.  They  seek  the  unknown  coefficient  through  minimization  of  the 
corresponding  “residual  fluxes”,  {<7,  —  7Vui}flj.  They  consider  minimizing  £,  where 


£  =  /  |<7,  —  7Vu,|2  dx. 

1=1 


The  idea  again  is  that  £  is  non-negative,  and  that  if  £  =  0,  then  the  conductivity  7  is 
consistent  with  the  given  boundary  measurements. 

For  the  numerical  implementation  of  this  algorithm  we  use  finite  element  approxima¬ 
tions  to  solve  the  boundary  value  problems  in  steps  (2. a)  and  (2.b).  We  take  piecewise 
linear  test  and  trial  functions  on  a  uniform  triangulation  of  the  domain  Q.  7  is  taken  to  be 
constant  on  each  element  T},  j  —  1,  •  ■  •  ,  N.  In  approximating  step  3  we  simply  minimize 
the  functional  F  with  fixed  {uj/hj  and  a ,  =  7V17,  1  <  i  <  K,  coming  from  the  finite 


1  We  use  the  least  squares  error  in  our  stopping  criterion. 


element  solutions  of  (2. a)  and  (2.b).  On  each  triangle  the  new  value  of  7  (in  the  scalar 
case)  is  simply  given  by 

(4  6)  7  -7 

Inew  -  laid  {  K  ) 

Di= 1 JT,  lVu«r  dx 

Since  u,  and  17  are  piecewise  linear,  the  integrals  are  not  necessary  and  we  get 

y>A'  |Vj;  ,2  1/2 

(4.7)  7 new  =  laid  ( 1  —  '  )  ,  on  each  triangle. 

E,=i  lv«.r 

As  a  measure  of  how  good  an  approximation  7 neu>  is  we  compute  the  LSE( least  squares 
error). 

LSE  =  y '[j  {7|Vu,|2  +  7-1|<7i|2}  dx  -  J  ds] 

(4.8) 

=  Y][f  {7neu>|Vu,|2  +  7~eWhoidVvi\2}  dx  -  f  <f>tip,  ds] 

,=  i  J a  Jr 

(remember  that  a,  =  7o/d  Viq). 

For  the  anisotropic  error  computation,  we  similarly  get 

(4.9)  LSE  =  y^[^{(7ncu,Vu,)rVui  +  ('f~'W~foidS7vi)T(joldVvi)}  dx  -  J  fcfc  ds]. 

5.  Two-Dimensional  Computation.  As  a  test  problem  we  consider  the  layered 
conductivity 


;’y)  =  {  0.5,  I 


2,  if  y  >  1/2 
0.5,  if  y  <  1/2 


on  the  d>  nain  Q  =  [0, 1]  x  [0, 1].  We  seek  to  reconstruct  this  7  from  three  sets  of  boundary 
measurements.  First  we  concoct  three  test  functions  U{,i  —  1,2,3  which  satisfy 

V  •  (7 Vii{)  =  0  in  Q 

(across  y  =  1/2,  this  requires  that  Ui  as  well  as  its  conormal  derivative  7%*-  be  continuous): 


“i(*,y)  =  j  J.  4 
l  4y, 

«2(^,y)  =  x 


«3 (x,y)  = 


y  +  1-5,  if  y  >  1/2 
4y,  if  y  <  1/2 

xy,  if  y  >  1/2 

4xy  —  1.5x,  if  y  <  1/2. 


. 


We  have  run  the  alternating  direction  algorithm  to  reconstruct  the  conductivity  in  both 
the  isotropic  and  the  anisotropic  cases  using  the  same  three  sets  of  boundary  measurements 

K>^)- 

In  the  actual  computations  to  solve  the  boundary  value  problems,  we  use  the  finite 
element  package  Modulef.  The  regular  finite  element  mesh  is  shown  in  figure  1.  Dirichlet 
and  Neumann  boundary  data  are  computed  from  the  exact  solutions  u,  and  prescribed 
pointwise  on  the  boundary  nodes  (for  Dirichlet  conditions)  and  as  an  average  value  on  the 
boundary  edges  (for  Neumann  conditions).  The  linear  systems  of  equations  are  solved  by 
Crout  decomposition. 

In  the  case  of  anisotropic  conductors  we  face  a  problem  of  nonuniqueness.  One  way 
to  avoid  the  nonuniqueness  issue  is  to  convert  the  matrix  conductivity  to  an  appropriate 
scalar.  In  the  graphic  display  shown  here  we  have  simply  taken  ^/det  7  as  an  isotropic 
approximation  at  the  final  stage,  but  we  recognize  that  this  is  an  ad  hoc  approach  which 
is  not  necessarily  best  possible.  For  the  graphics  display,  we  use  a  grid  point  in  each 
triangle,  not  the  nodes  of  the  original  triangulation.  (Figure  2  shows  the  grid  points  for 
the  graphics.) 

Figures  3  and  4  show  the  results  of  our  computation.  As  seen  both  the  isotropic  and  the 
anisotropic  reconstructions  detect  some  jump  near  the  line  y  =  j,  though  the  isotropic  re¬ 
construction  detects  a  sharper  jump  than  the  anisotropic  one.  Most  noticably  the  isotropic 
reconstruction  appears  to  become  oscillatory  after  a  certain  number  of  iterations,  whereas 
the  anisotropic  reconstruction  displays  no  such  oscillations  and  seems  to  improve  with 
increasing  number  of  iterations.  This  is  not  totally  unexpected,  since  the  reconstruction 
algorithm  using  anisotropic  7  was  derived  from  the  isotropic  algorithm  through  the  elimi¬ 
nation  of  highly  oscillatory  7  by  relaxation  (or  (7-convergence).  See  Ivohn  and  Vogelius  [8]. 
In  both  the  isotropic  and  the  anisotropic  cases,  the  least  squares  error  decreases  rapidly. 
As  seen  the  anisotropic  reconstruction  brings  about  smaller  least  squares  error  than  the 
isotropic  one  when  we  compare  the  results  at  the  same  iterate.  (See  table  1.) 

To  examine  how  close  the  reconstruced  conductivity,  yc,  is  to  the  exact  conductivity, 
7e,  we  have  computed  |j 7°  -  7e||  using  both  the  Lj  and  L2  norms.  (See  figure  5  and  table 
2.)  In  this  particular  example,  we  see  that  the  anisotropic  L\  error  decreases  and  stays 
almost  unchanged  after  about  200  iterations,  while  the  isotropic  L\  error  rapidly  decreases 
at  the  very  early  sfage  of  computation  but  starts  increasing  when  7  becomes  oscillatory 
(at  about  the  50' h  iterate)  and  remains  almost  unchanged  after  about  200  iterations.  The 
two  error  curves  intersect  at  about  the  120<ft  iterate,  after  which  the  isotropic  Lx  error 
curve  lies  significantly  above  the  anisotropic  one.  We  also  have  computed  the  Z2  error 
and  the  results  are  qualitatively  similar  to  those  obtained  using  the  L\  norm,  although  the 
crossing  happens  at  a  later  stage  and  is  less  dramatic. 

6.  One-Dimensional  Computation.  We  now  incorporate  into  the  algorithm  the 
fact  that  7  is  a  function  of  x2  only.  We  take  periodic  boundary  conditions  in  the  X] 


a 

S5SSS1S8S 

Si; 

direction  and  for  simplicity  we  change  the  domain  S2  to  be  [0,  27rj  x  [—1, 1J.  Assume  the 
solution  u  is  represented  by  a  sine  Fourier  series 

OO 

u(x i ,  X2)  =  Y^  at(i2)  sin  kx  1 

k=l 

(no  cosine  terms).  To  minimize  fn  7|Vu[2  dx  is  equivalent  to  find  the  solutions  of 

~(l<x'k)'  +  =  0,  k  =  1, 2, • ■ • 

with  appropriate  boundary  data  and  to  minimize  /n  7-1|<t|2  dx  subject  to  V  •  a  =  0  is 
equivalent  to  find  the  solutions  of 

-h~lk~2a'k)'  -f  7_1crit  =  0,  =  1,2,-- 

with  appropriate  boundary  data.  The  latter  follows  since  a  is  represented  as 

/  Y2'  \ _ \ 


<r(x\ ,  X2  )  — 


Y  -Ta'k(x2)cOskx1 


=  *=1 
00 


Y  ak(x2) sin  kxl 


We  seek  to  minimize  the  truncated  funtional 

K  fl 


(6.1)  F(7,{afc}^Li,{«7fe}^=1)  =  Y  /  h(k2\<*k]2  +  K-]2)  +  7  \k  2Wk ? +  [vk?))dx 

k=i 

over  all  a*  and  a  *  with  given  boundary  data  <>*(±1)  and  <rt.(±i). 

The  Alternating  Direction  Algorithm  for  the  one-dimensional  computation  becomes: 

1.  Pick  an  initial  guess  for  7. 

2.  For  £  =  !,•••  ,K 


2. a  Solve 


2.b  Solve 


-  ha'k)'  +  k2l<*k  =  0  in  (-1,1) 
a*(  —  1)  and  or*(l)  given. 

-  (7~1fc-2o-Jt)' +  7_1^fc  =  0  in  (-1,1) 
(7*..(-l)  and  <7jt(l)  given. 


3.  Update  7  by  minimizing  F(~/,  {  c*  *}*._! ,  {crfc}£=1)  for  fixed  {a*}*'-!  and  {cjt}t=1- 

4.  If  the  reconstruction  is  satisfactory,  exit;  else  go  to  step  2. 

For  the  numerical  implementation  of  this  one-dimensional  algorithm  we  use  piecewise 
linear  finite  elements  on  a  uniform  mesh  to  solve  the  boundary  value  problems  (2. a)  and 
( 2. b ) .  We  implement  the  method  as  described  above  for  two  component  materials.  We 
also  implement  the  anisotropic  analogue  of  this  algorithm  based  on  volume  fractions. 


WWW w? 


mmrnm* 


6.1  Material  composed  of  two  components.  For  the  minimization  of  F  we  con¬ 
sider  conductivities  7  that  correspond  to  a  two  component  material,  i.e.,  7  takes  one  of 
two  values  7^  or  7^  on  each  element  I,,  i  —  1,  •  •  •  ,  N.  We  seek  to  minimize 


Zb  f  +  Y^[Q'k}2)dx  +  ~iT1  f  (Zfc_2bt]2  +  Z^2)dx]’ 

1=1  Jl‘  Jt-l  fc=l  k=l  k= 1 


which  amounts  to  minimizing 


K 

hY^Jji^iak)2  4-  \*X)dx  +  7,-1  E/(*-2K]2  +  [ok}2)dx 


on  each  element 

We  compute  (6.2)  with  7,  replaced  by  7^  and  7^  respectively,  compare  the  resulting 
values  and  select  as  the  new  conductivity  the  one  which  gives  rise  to  the  smaller  energy. 
Step  3  may  now  be  written: 


3.  For  i  =  1,  •  •  •  ,  N 


K  K 

di=7(1)Z  /  (*2[^]2  +  K.]2)<fx  +  7(1)  *Z  f  (k~2[a'k]2  +  [ak]2)dx 

\Jli  k=iJli 

(6.3)  rf2  =  7(2)Z  f  (A:2  [a*]2  +  [a'k]2)dx  +  7(2)  *  /  (k~2  [a'k]2  +  [a k]2 )  dx 

Jr, 


If  d\  <  </2,  then  7*  =  7*^; 

else  7 i  =  7^. 


The  least  squares  error  is  computed  by  summing  all  d\  or  cf2,  whichever  is  smaller  on  each 
element,  and  finally  subtracting  the  boundary  integral. 

For 

f  2,  for  —  1  <  x  <  0 
7(x)  =  < 

l  0.5,  for  0  <  x  <  1, 


c**(x)  and  <7 k(x),  1  <  k  <  K ,  can  be  computed  analytically 


ekx  +  e-kx 


ak(x^  tu  x 


k(ek  +  e  *) 


crk(x)  = 


ekx  _  e-kx 


2  —r - j-r,  for  -1  <  X  <  0 

(ek  +  e~k)  ~  ~ 


ekx  _  e-kx 


0.5  7— r - j—,  for  0  <  x  <  1. 

( ek  +  e~k )' 
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We  feed  the  algorithm  with  the  exact  boundary  data  ajt(±l)  and  crjt(±l),  1  <  k  <  K, 


ajt(  — 1)  =  ojt(l)  =  1  /k 


— 1)  =  - 


—ek  +  e  k 
ek  -f  e~k 


crjt(l)  =  0.5 


ek  +  e  k 


In  the  actual  computation  we  take  I\  =  10  and  N  =  200.  We  have  run  the  alternating 
direction  algorithm  on  this  problem  to  reconstruct  7  for  various  initial  choices.  In  each 
case  a  new  conductivity  was  obtained  after  one  iteration  and  then  it  remained  the  same 
in  all  subsequent  iterations.  The  least  squares  error  changed  from  the  first  to  the  second 
iterate,  since  the  computation  of  the  least  squares  error  involves  7 as  well  as  7 new . 

Here  are  some  results  with  different  initial  guess. 

1.  Initial  7  =  2.0  everywhere:  the  algorithm  picks  up  a  jump  at  x  —  0.29  with  least 
squares  error  0.01638.  Computed  conductivity  7C  is:  7C  =  2  for  x  <  0.29,  yc  =  0.5 
for  x  >  0.29. 

2.  Initial  7  =  0.5  everywhere:  the  algorithm  picks  up  a  jump  at  x  =  0.29  with  least 

squares  error  0.01638.  7C  =  2  for  x  <  0.29,  =  0.5  for  x  >  0.29. 

f  0.5,  for  -1  <  x  <  0 

3.  Initial  'y(x)  =  <  :  the  algorithm  picks  up  a  jump  at  x  =  0.58 

l  2,  for  0  <  x  <  1 

with  least  squares  error  0.08429.  -yc  =  2  for  x  <  0.58,  7C  =  0.5  for  x  >  0.58. 

4.  Initial  7  =  exact  7:  the  algorithm  preserves  the  initial  guess  (jump  exactly  at 
i  =  0)  with  least  squares  error  0.00144,  the  descretization  error  in  the  numerical 
solutions  of  the  boundary  value  problems. 

6.2  Volume  fraction  method.  The  method  discussed  in  the  previous  section  itself 
failed  to  detect  the  exact  location  of  the  jump  discontinuity.  The  method  would  be  equally 
unsuited  to  identify  any  oscillations  in  7.  We  now  consider  a  method  for  a  two  component 
material  based  on  the  concept  of  volume  fractions  (the  anisotropic  analogue  of  the  previous 
method). 

With  a  possibly  highly  oscillatory  layered  coefficient,  it  is  well  known  from  the  theory  of 
homogenization  that  the  conductivity  across  the  layers  approaches  the  harmonic  average 
c  =  (#7^)  +  (1  —  8)  7^  )  1 ,  whereas  the  conductivity  along  the  layers  approaches 

the  algebraic  average  m  =  07^  +  (1  —  6)^2\  here  8  is  the  so  called  volume  fraction  of 
material  7^,  i.e.,  the  infinitesimal  proportion  of  material  7^.  See  Bensoussan,  Lions, 
and  Papanicolaou  [2].  Clearly,  0  <  0(x)  <  1  and  we  note  that  if  either  8  =  0  or  1,  then 
c  =  m  =  7(2^  or  7^,  respectively. 

We  seek  to  minimize  the  truncated  functional 


(6.4)  F(0,  {«*}£=„ 


-f  mA:2[ajt]2  +  m  lk  2[<rJ.]2  +  c  1  [oq]2  }  dx. 


We  take  c  and  m  to  be  constant  in  each  element  Ii. 

Alternating  Direction  Algorithm  based  on  volume  fraction  9: 

1.  Pick  an  initial  guess  for  0. 

2.  For  k  =  !,•••,  I\ 


Solve 


2.b  Solve 


—  ( ca'k )'  +  k2motk  =  0  in  (  —  1, 1) 
ajt(  — 1)  and  ajt(l)  given. 

—  ( k~ 2  m~la'k)'  +  c~1<T/ 1  =  0  in  (  —  1, 1) 
crjt(-l)  and  au(l)  given. 


3.  Update  9  by  minimizing  F(9,  {at}£Lj,  {<7i}*=1)  for  fixed  and  {<7jt}£_1. 

4.  If  the  reconstruction  is  satisfactory,  exit;  else  go  to  step  2. 

In  updating  9  in  step  (3),  we  minimize  the  nonlinear  functional 

(6.5)  Ft{c)  =  Y(c  f  [a'k]2  dx  +  mk2  f  [ak]2  dx  +  m~1  k~2  f  [cr'k]2  dx  +  c~*  [  [ck]2dx), 
Jit  Jli  Ju  Jli 

with  respect  to  c  on  each  element  Jj,  and  we  compute  9  from  the  formula 

9  =  (c-1  —  7^  1)/(7 1  —  ). 

Note  that  F,  is  a  functional  of  c  only,  since  m  can  be  expressed  in  terms  of  c 

(U+  (2)  7(1)7(2) 

c 

In  our  computation  to  minimize  Fj,  we  use  Newton’s  method  together  with  the  Golden 
search  rule.  We  take  two  Newton  iterates.  If  the  second  iterate  C2  falls  between  7U)  and 
7<2) ,  then  we  compare  F, ( C2 ) ,  Fi( 7'^),  and  F,(~f(-2'1 )  and  choose  the  argument  corresponding 
to  the  smallest  value  as  the  updated  value  for  c.  If  the  second  iterate  C2  falls  outside  the 
range,  then  we  use  the  Golden  search  algorithm  (ten  steps)  to  compute  an  updated  value 
for  c. 

The  least  squares  error  is  computed  by  summing  the  F,’s  for  all  i  and  finally  subtracting 
the  boundary  integral.  Here  we  use  both  0oid  and  9new  in  the  sense  that  9ntw  is  used  to 
compute  F,  and  90m  is  used  to  solve  the  boundary  value  problems.  For  our  test  problem 
we  take  the  same  exact  solutions  as  in  the  previous  section,  i.e.,  7^  =  0.5  and  7^  =  2, 
and  in  terms  of  volume  fraction 


■It 


0,  for  —1  <  x  <  0 
1,  for  0  <  x  <  1. 


We  have  run  the  algorithm  on  this  problem  to  reconstruct  6  (again  with  K  =  10, 
N  =  200  and  exact  boundary  data)  with  various  choices  of  initial  guess;  the  results  however 
are  almost  independent  of  initial  guess.  This  method  reconstructs  9  fairly  well  succeeding 
to  detect  the  jump  discontinuity.  Figure  6  shows  the  results  obtained  from  this  method 
with  initial  guess  9  =  0.5.  The  least  squares  error  decreases  very  rapidly  at  the  very 
early  stages  of  the  computation.  We  have  also  computed  ||#c  —  0e||z,i,  where  9C  is  the 
reconstructed  volume  fraction  and  9e  is  the  exact  volume  fraction.  (See  table  3.)  As  seen 
the  L 1  error  decreases  as  the  number  of  iterations  increases. 

We  have  tested  the  same  problem  on  meshes  of  various  dimensions.  Figure  7  shows  the 
reconstruction  results  when  N  =  25.  We  see  that  the  least  squares  error  is  much  smaller 
when  N  is  larger.  For  example,  the  least  squares  error  when  N  =  25  is  0.07249  at  the 
50t,>  iterate  while  it  is  0.00122  when  N  =  200  at  the  same  iterate.  Surprisingly  enough, 
the  L\  error  is  smaller  when  N  is  smaller.  For  example,  the  L\  error  when  N  =  25  is 
0.05641  at  the  50f/l  iterate  while  it  is  0.0.08921  when  N  =  200  at  the  same  iterate.  This 
is  a  phenomenon  which  we  have  come  across  in  several  of  the  computations  -  basically  it 
asserts  that  there  is  nothing  to  be  gained  (occasionally  something  to  be  lost)  by  using  a 
too  fine  mesh  to  identify  a  simple  discontinuous  coefficient.  See  figure  8  and  table  3  for 
the  comparison  of  the  least  aquares  errors  and  the  L\  errors. 

6.3  Perturbed  data.  Let  us  consider  a  problem  to  reconstruct  9 ,  slightly  perturbed 
from  the  one  in  the  previous  section,  i.e.,  a  problem  which  corresponds  to  the  exact  9 

(  —  jx  —  for  -1  <  x  <  —0.5 
0,  for  -0.5  <  x  <  0 

9(x)  =  < 

1,  for  0  <  x  <  0.5 

k  —  jx  +  for  0.5  <  x  <  1. 

We  do  not  compute  the  boundary  data  analytically  for  this  problem.  Instead  we  set 
o*;(±l)  =  1  fk  for  k  =  1,  •  •  •  ,  K  and  solve  the  boundary  value  problems 

—  (coIl)'  +  mfc2ajt  =  0  in  (  —  1,1) 

(6.6)  *  ' 

ajfc(-l)  =  «*(  1)  =  1/fc 

numerically,  using  piecewise  linear  finite  elements  on  a  mesh  with  N  =  800.  As  data  for 
(7*.-,  we  use 

<t*(-1)  =  c(— l)a'M(  — 1) 
a*(l)  =  c(l)a'M(l), 

where  a'k  h(x)  is  the  finite  element  approximation  of  a'k(x). 

We  have  run  the  Alternating  Direction  Algorithm  based  on  volume  fractions  on  this 
problem.  The  reconstruction  results  are  shown  in  figure  9  and  table  4  with  initial  guess 


9  =  0.5.  Almost  the  same  conclusions  can  be  drawn  from  this  computation  as  in  the  case 
of  the  previous  problem  except  that  the  L1  error  is  larger  for  smaller  N  in  this  case. 

Here  are  some  observations: 

1.  The  algorithm  reconstructs  8  fairly  well  with  reasonably  small  least  squares  error. 

2.  The  results  are  almost  independent  of  initial  guess  for  9. 

3.  The  least  squares  error  decreases  rapidly  as  N  gets  larger  when  compared  at  the 
same  iterate. 

4.  The  L\  error  in  the  computed  volume  fraction  decreases  as  the  number  of  iterations 
increases. 

For  comparison,  we  tried  the  algorithm  in  section  6.1  on  this  problem  to  see  whether 
it  would  detect  any  oscillations.  It  picked  up  a  single  jump  after  one  iteration  and  the 
solution  remained  the  same  in  all  subsequent  iterations.  The  location  of  the  computed  jump 
depends  on  the  initial  choice  of  7,  however  the  method  failed  to  pick  up  any  oscillations. 
The  location  of  jump  is  slightly  different  from  the  one  in  section  6.1  with  the  same  initial 
guess;  the  least  squares  error  is  naturally  much  bigger.  For  example,  with  the  initial  guess 
7  =  0.5  everywhere,  the  algorithm  picked  up  a  jump  at  x  =  0.25,  i.e.,  7C  =  2  for  x  <  0.25, 
7C  =  0.5  for  x  >  0.25  with  least  squares  error  0.43973. 

6.4  Piecewise  quadratic  finite  elements.  As  seen  we  have  obtained  quite  satis¬ 
factory  reconstructions  of  0  in  the  previous  two  sections.  To  analyze  the  effect  of  higher 
accuracy  we  also  tried  to  use  piecewise  quadratic  finite  elements  in  solving  the  boundary 
value  problems  (2. a)  and  (2.b).  The  coefficients  c  and  m  are  treated  as  piecewise  linear  (not 
necessarily  continuous).  In  the  step  to  update  9  it  might  seem  natural  to  minimize  over 
all  piecewise  linear  9.  This  however  produces  a  very  complicated  updating  step.  Instead 
we  have  decided  to  use  one  of  the  following  two  options  :  (1)  keep  9  piecewise  constant 
as  before,  or  (2)  minimize  the  functional  with  respect  to  8  pointwise  at  the  nodes  of  each 
element  and  let  9  elementwise  be  the  linear  interpolation  (not  necessarily  continuous). 

We  have  run  the  Alternating  Direction  Algorithm  based  on  the  volume  fractions  on 
the  problems  of  section  6.2.  Figure  10  and  table  5  show  the  results  of  this  computation.  In 
our  first  computations  we  reconstructed  the  volume  fractions  using  option  (1).  The  results 
are  almost  the  same  as  with  piecewise  linear  elements  in  the  sense  of  L\  error,  of  course 
the  least  squares  error  is  significantly  reduced.  In  the  actual  computation  of  the  problem, 
we  got  a  negative  least  squares  error  in  the  computation  with  N  =  200;  this  is  due  to 
round-off  error.  We  computed  the  same  problem  in  double  precision,  which  brought  about 
a  very  small,  but  positive  least  squares  error. 

We  have  also  used  option  (2)  on  the  same  problem.  The  results  look  reasonable  ex¬ 
cept  near  the  jump  discontinuity.  This  is  not  unexpected  since  the  minimization  of  the 
functional  requires  pointwise  approximations  to  the  derivatives  of  the  solutions  to  the 
boundary  value  problems.  These  approximations  are  clearly  not  very  good  near  the  jump 
discontinuity. 


Finally  we  computed  with  the  perturbed  data  from  section  6.3.  We  used  option  (li¬ 
the  results  are  shown  in  figure  10  and  table  6. 

In  summary  it  seems  that  the  increased  accuracy  in  the  finite  element  solutions  has 
very  little  effect  on  the  L\  error  of  the  computed  solutions.  Based  on  this  observation  and 
considerations  regarding  simplicity,  it  seems  entirely  reasonable  to  solve  the  finite  element 
problems  using  only  piecewise  linears.  We  do  not  know  whether  this  conclusion  is  also 
valid  for  more  smoothly  varying  9.  In  that  case  one  might  suspect  that  higher  accuracy 
could  improve  the  L\  error. 
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Table  1 


of  iterations 

least  squares  error 

least  squares  error 

(isotropic) 

(anisotropic) 

1 

0.43236 

0.26953 

5 

0.C0781 

0.00181 

10 

0.00496 

0.00081 

20 

0.00274 

0.00056 

30 

0.00161 

0.00044 

40 

0.00109 

0.00037 

50 

0.00081 

0.00033 

100 

0.00033 

0.00024 

150 

0.00023 

0.00019 

200 

0.00022 

0.00015 

250 

0.00019 

0.00011 

300 

0.00019 

0.00010 

350 

0.00019 

0.00010 
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Table  2 


number  of  iterations 

Ll  error 

Ll  error 

L2  error 

L2  error 

(isotropic) 

(anisotropic) 

(isotropic) 

(anisotropic) 

5 

0.230062 

0.247373 

0.3588591 

0.366115 

20 

0.178497 

0.233846 

0.295077 

0.364288 

50 

0.166697 

0.215210 

0.283006 

0.360893 

100 

0,185301 

0.197864 

0.307787 

0.352130 

150 

0.200615 

0.186805 

0.326579 

0.346145 

200 

0.208475 

0.182351 

0.336369 

0.343204 

250 

0.212179 

0.180737 

0.341201 

0.342213 

300 

0.213886 

0.181449 

0.343713 

0.342388 

350 

0.214690 

0.183172 

0.345156 

0.34012 

Reconstructed  volume  fraction 
50-th  Iterate  when  N-200 
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Table  3 


number  of  iterations 
Figure  8 


of  iterations 

least  squares  error 

least  squares  error 

Li  error 

Litnor 

(N=2S) 

(N=200) 

(N=25) 

(N=200) 

10 

0.07325 

0.00199 

0.09785 

0.13418 

20 

0.07271 

0.00144 

0.07809 

0.11360 

30 

0.07257 

0.00131 

0.08726 

0.10274 

40 

0.07252 

0.00125 

0.06088 

0.09508 

50 

0.07249 

0.00122 

0.05641 

0.08921 

Reconstructed  volume  fraction 
50-th  Iterate  when  N-25 


Reconetructed  volume  fraction 
50-th  Iterate  when  N-200 
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of  iterations 

least  squares  error 

least  squares  error 

Li  error 

Tjerror 

(N=25) 

( N  =  2  00 ) 

( N  =  2  5 ) 

(N=200) 

10 

0.08502 

0.00246 

0.25482 

0.21760 

20 

0.08439 

0.00169 

0.22025 

0.I9I48 

30 

0.08117 

0.00156 

0.19988 

0.17979 

40 

0.08403 

0.001  18 

0.18331 

0.171 10 

50 

0.08392 

0.00144 

0.16895 

0.16420 
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Table  5:  Piecewise  quadratic  finite  element  results  on  problem  in  sec’  ■  6.2 


When  option  (1)  is  used: 


number  of  iterations 

least  squares  error 

least  squares  error 

L i error 

Zierror 

(14=25) 

(N=200) 

(N=25) 

(N=200) 

(single  precision) 

(double  precision) 

(single  precision) 

(double  precision) 

10 

0.00116 

0.00091 

0.11987 

0.13418 

20 

0.00067 

0.00031 

0.10005 

0.11362 

30 

0.00055 

0.00017 

0.08851 

0.10275 

40 

0.00050 

0.00011 

0.07993 

0.09509 

50 

0.00047 

0.00008 

0.07331 

0.08921 

When  option  (2)  is  used 

number  of  iterations 

least  squares  error 

least  squares  error 

terror 

Z^error 

(N=25) 

(N=200) 

(N=25) 

(N=200) 

(single  precision) 

(double  precision) 

(smi-le  precision) 

(double  precision) 

10 

0.00344 

0.00105 

0.11689 

0.13454 

20 

0.00276 

0.00043 

0.09963 

0.11525 

30 

0.00213 

0.00027 

0.09254 

0.10544 

40 

0.00167 

0.00020 

0.08457 

0.09858 

50 

0.00155 

0.00015 

0.07799 

0.09327 

Table  6:  Piecewise  quadratic  finite  element  result.'  on  problem  in  section  6.3 


number  of  iterations 

least  squares  error 

least  squares  error 

Li  error 

L\  error 

(N=25) 

(N=200) 

(N=25) 

(N=200) 

(single  precision) 

(double  precision) 

(single  precision) 

(double  precision) 

10 

0.00326 

0.00129 

1.20264 

0.21753 

20 

0.00174 

0.00039 

0.17177 

0.19132 

30 

0.00162 

0.00025 

0.15869 

0.17944 

40 

0.00144 

0.00019 

0.14894 

0.17061 

50 

0.00140 

0.00012 

0.14094 

0.16356 

